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Abstract The Sudakov veto algorithm for generating emission and no-emission 
probabilities in parton showers is revisited and some reweighting techniques are 
suggested to improve statistics by oversampling in specific cases. 



1 Introduction 

Sudakov [T| form factors or no-emission probabilities, are used in all parton shower 
(PS) programs, to ensure exclusive final states and, at the same time, to resum 
leading (and sub-leading) logarithmic virtual contributions to all orders. The way 
these enter into the shower is through the ordering of emissions, typically in trans- 
verse momentum or angle. In the standard case we have an inclusive splitting 
probability of a parton % into partons j and k given by 
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dVi,jk(t,z) = -lPijk(z)dtdz, (1) 

Z7T 

where t = logpj_/ Aq CD i s the ordering variable, z represents the energy sharing 
between j and k, and where we have integrated over azimuth angle in the Altarelli- 
Parisi splitting function Pijk(z). Starting from some maximum scale to we then 
want to find the exclusive probability of the first emission, which we get from the 
inclusive splitting probability by multiplying with the probability that there is no 
emission before the first emission, 

dPij* (t. z) = ^P iljk (z)dtdz x A(t ,t). (2) 

Here, A(to,t) is this no-emission probability, or the Sudakov form factor, given by 

A(t , t) = exp {- J ° dt'dz^Pijkizfj . (3) 

In principle the Sudakov form factor can be calculated analytically. However, 
often the integration region in the z-integral can be non-trivial, and most PS 
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programs today prefer to calculate it numerically using the so-called Sudakov veto 
algorithm [2]. The trick here is to find a simple function which is everywhere 
larger than Pijk and which is easy to integrate, and by systematically vetoing 
emissions generated according to this overestimated function, a correct no-emission 
probability is obtained. 

The Sudakov veto algorithm (SVA) is normally used for purely probabilistic 
processes, but recently it has been generalized to also be used in cases where the 
function being exponentiated is not positive definite OH]. 

In this article we shall investigate other modifications of the Sudakov veto 
algorithm, where we try to increase the statistical precision in some special cases 
by oversampling techniques, but we will also briefly discuss the issue of negative 
contributions to splitting functions. 

In section [3] we will investigate the usage of the SVA in CKKW-like 016] 
algorithms, where parton showers are matched with fixed order inclusive matrix 
elements (MEs). Here, the SVA is used to make the inclusive MEs exclusive by 
multiply them with no-emission probabilities taken from a parton shower. One 
problem with this procedure is that every event prduced with the matrix element 
generator is either given the weight zero or one, which becomes very inefficient if 
the cutoff used in the ME-generation is small. We will find that by introducing 
oversampling, a weight can be calculated which is never zero, but nevertheless 
will give the correct no-emission probability. In section [3] we will also discuss an 
extension of the CKKW algorithm to include next-to-leading order (NLO) MEs [7] 
where the SVA is used to extract fixed orders of a s form the parton shower to avoid 
double counting of corresponding terms in the NLO calculation. 

Then, in section [5] we will consider cases where a parton shower includes dif- 
ferent competing processes, where some of them are very unlikely. This is the case 
in e.g. the PYTHIA parton shower, where photon emissions off quarks are included 
together with standard QCD splittings. Since «em is much smaller than a s it is 
very time consuming to produce enough events containing hard photons to get 
reasonable statistics. We shall see that a naive oversampling of the photon emis- 
sions has unwanted effects on the total no-emission probability, and that a slightly 
more involved procedure is needed. The method presented is different from the 
one introduced by Hoche et al. in [8], but is equally valid. It turns out that both 
these methods can be used to include negative terms in the splitting functions. 

But first we shall revisit the derivation of the SVA, as we will use many of the 
steps from there when we investigate the different oversampling techniques. 

2 The Sudakov Veto Algorithm 

Here we follow the derivation found in [2] and [9]. Although we normally have 
competing processes, we will first simplify the notation by just considering one 
possible splitting function, which we will integrate over z, 





The no-emission probability simply becomes 



A(t ,t c ) = e 




(5) 
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If r(t) can be integrated analytically, and if the primitive function, f has a 
simple inverse, it is easy to see that we can generate the t-value of the first emission 
by simply generate a random number, R, between zero and unity and obtain 

f ° dt'r(t')A(t ,t') = R [ ° dt'r(t')A(t ,t') 
Jt Jo 



l-A(t o ,t) = R(l-A(t o ,0)) 
1-exp ( f(t) — r(t )) = R 
t = f~ 1 (f(t ) + \og(l-R)), 



(6) 



where we have assumed that r(t) is divergent for small t, such that A(t, 0) = 0, 
an assumption we will come back to below. 

Now, in most cases the integration of r is not possible to do analytically, and 
if it is, the inverse function may be non-trivial. This is the case which is solved by 
the SVA. All we need to do is to find a nicer function, _T, with an analytic primitive 
function, which in turn has simple inverse, such that it everywhere overestimates 

r, 

f(t)>r(t), vt. (7) 

With this function we can construct a new no-emission probability 



A(to,t c ) 



■ C p ^ dt 



(8) 



which is everywhere an underestimate of A(to,t c ), and we can generate the first t 
according to it. As in the standard accept-reject method, we now accept the gen- 
erated value with a probability r(t)/r(t) < 1. However, contrary to the standard 
method, if we reject the emission, we replace to with the rejected i-value before we 
generate a new t. Loosely speaking, we have underestimated the probability that 
the emission was not made above t, so we need not consider that region again. We 
now continue generating downwards in t until we either accept a lvalue, or until 
the generated t drops below t c at which point we give up and conclude that there 
was no emission above t c . 

To see how this works more precisely, we look at the total probability of not 
having an emission above t c , 



Via 



5>„, 



0) 



where V n is the probability that we have rejected n intermediate ^-values. To start 
with, we have 



To = A{t ,t c ) 



dtf(t)A(t ,t) 



r(t) 



A(to,t c ) [ ° dt[f(t)-r(t)] 



A(t,t c ) 



(10) 
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where for Vi we first have the probability that we generate a value t and then 
throw it away with probability 1 — r(t)/r(t) and then the probability that we do 
not generate anything below t. Similarly we get 



V 2 



rto 



dtxr^AitoM) 



1 



nti) 



nti) 



dt 2 r(t 2 )A(t 1 ,t2) 



EM 



A{t 2 ,t c ) 



(ii) 



from noting that A(to, t\)A{t\, t 2 )A(t 2 , t c ) = A(to,t c ) and that the nested integral 
can be easily factorized. 

This is easily generalized to an arbitrary number of vetoed emissions and we 

get 

oo oo , to s n 

J2 T - = J2 A (M^( dt[f(t)-r(t)]\ 



n=0 

A(t ,t c )ei: 
A(t ,t c ), 



(12) 



which is the no-emission probability we want. 

Here we have ignored the additional variables involved in the emissions, ft is 
easy to see that if we have a simple overestimate of the splitting function 



P(t,z)>^P(z), Vt,z 



we can construct our r as 



f(t) = [ 

J z 



*„„(t) 



P(t,z)dz 



(13) 



(14) 



2„in(t) 



with z max (i) > z m ax(i) and i m in(i) < Zmin(i) overestimates the integration re- 
gion in z. For each t we generate, we then also generate a z in the interval 
[zmin(i), ^max(i)] according to the probability distribution 



V{z) 



r(t) 



and we then only keep the emission with the probability 

^P(z)Q(z - z min (<)e( W(*) - z)) 



P(t,z) 



(15) 



(16) 



Although the formulae become more cluttered, it is straight-forward to show, by 
going through the steps above that this will give the correct distributions of emis- 
sions. 
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If we now go back to eq. ([6j> , we there assumed that r diverges at zero such 
that A(t,0) = 0. This is, of course, not necessarily the CcLSG, cLS pointed out in [3] . 
However, here we will only concern ourselves with emissions above some cutoff 
t c > 0, and we can always add to our overestimate a term which is zero above t c 
but which diverges at t = 0. The fact that nowhere in the veto algorithm does 
anything thing depend on the form of P below t c , means that we do not even need 
to specify how it diverges, it is enough to assume that it does. 

In parton showers one always have many different emission probabilities. We 
typically have several possible emissions for each parton. The SVA is easily ex- 
tended to this situation by noting that the no-emission probability factorizes for 
the sum of different splittings, 

A(t ,t c ) = e~-C £» r " (t)dt = Y[A a (to,t c ). (17) 

a 

The first emission is then generated with a t-value given according to 

dV = ^2r a (t)A(t ,t)dt, (18) 

a 

and then randomly selecting one of the processes with weight 

It is easily shown that generating one t for each possible splitting according to 

dVa = r a (t a )A a (t ,t a )dta, (20) 
and then selecting the splitting which gave the largest t a , gives the same result. 



3 Reweighting in CKKW-like procedures 

In CKKW we want to multiply partonic states generated by a ME generator with 
Sudakov form factors. We take the partonic state and project it onto a parton- 
shower history. An n-parton state is then reconstructed as a series of parton shower 
emissions with emission scales {ti,...,t n } and the corresponding intermediate 
states {So, ■ ■ ■ , S n } where S n is the one generated by the ME. We then want to 
multiply by the no-emission factors 

- f !i dtr-At) 

Ai(ti,t i+1 ) = e Jt <+i (21) 

where Jj is the sum of the splitting functions from the partons in state Si- 

What we can do is to simply put the state Si into the parton shower program 
and ask it to generate one emission starting from the scale U. If the generated 
emission has a scale t > tj+i we throw the whole partonic event away and ask the 
ME generator to produce a new state. The probability for this not to happen is 
exactly Ai(ti, ij+i) and the procedure corresponds to reweighting the ME state 
with the desired no-emission probability. 
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The problem is that if the ME state corresponds to very low scales, we will 
throw away very many events, which is very inefficient and may result in poor 
statistics. 

A way to improve this situation is to introduce a boost factor for the splittings 
■Ti — > Pi = CTi, with C > 1, and multiply the overestimate Pi with the same 
factor. As before this just gives a simple overestimate of the splitting function, 
which we know how to handle from section[2] But rather than throwing an emission 
away with an extra probability 1/C (and not veto the event) we can always reject 
the emission but multiply the whole event with a weight 1 — 1/C. The total weight 
of the event will then be the sum of all possible ways we can veto a generated 
emission (we here assume that the normal rejection procedure has already been 
applied) 

( Wi ) = jT(l - l/C) n Ai(ti,t i+ i)^ ( I 1 dtPi(t)) 



00 if r u 

71 = \ J U + 1 



dt [r z (t) - r z (t)] J 

= Ai(ti,t i+ i) (22) 

In this way we get the right weight but we never throw away an event. 

In the NLO version of the CKKW-L [7J we also want to calculate the integral of 
the splitting function, J f ' dtr(t)i which is used as a way of subtracting the fixed 
first orders result from the exponentiation and then replace it with the correct 
NLO result. The way this was done in [7J was similar to the procedure above. The 
shower is started, and each emission is vetoed, but the number of emissions above 
ti+\ was counted and it was noted that the average number of vetoed emissions is 
given by 

(n) = V nAi{ti,U +1 )—. ( / dtr(t)i) 

= { ' dtTi(t) x ]T Ai(ti,ti+i) * ( f 1 dtr^t) \ 

= /" 



l-ti 

dtr^t) (23) 



To factor out other fixed order terms we note also that 

(n(n-l)-..(n-m)> = —. I / dtTtft) I , (24) 




which means we can pick out higher-order terms in a a used in the shower. 

Again, the statistics can become a bit poor if most events yield the weight 
zero (which is the case if for large merging scales when the no-emission probability 
is close to unity), and only a few have non-zero values. We can instead again 
introduce the boost factor, C, and rather than simply counting the number of 
emission we take the weight n/C. 
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We note that in general this C need not be a simple constant, it can be a 
function of the scale (or any other variable in the splitting.) This is used in the 
NLO version of CKKW-L, where the leading order a s term in the expansion of 
the no-emission probability is needed at a fixed renormalization scale, hr, while in 
the shower we have a coupling running with the transverse momentum. Therefore, 
rather than counting the the number of emissions, we sum up ratios of fixed and 
running a s for the emissions which are generated and discarded. Introducing 

= / P(z)dz = r(t) (25 

i 2mill (() 2n «s(t) 

We then get on the average a weight 



(w) = V A(t , to) / dtir(ti) ■■■ dt n r(t n ) V 

n=0 Jt c Jt Vi=i as[ - U) 



to 

&T H (*), (26) 



by working a bit on the symmetrized nested integrals of different functions and 
we get what we desired. To obtain higher powers of the integral with fixed q s , it 
is easy to show that e.g. the average sum of triplets 

Qs(Mfl) Qs(mh) Qs(Mfl) , 2r 
.-f-' a s {U) a s (tj) a s {t k ) 

will give 



1_ 
3! 



to \ 3 

dtr R (t)\ . (28) 



We also note that for initial-state splittings, the integral over splitting functions 
also contain ratios of parton density functions, 

fb(x/z,t) 

Ux,t) ' (M) 

where a is the incoming parton before, and b is the parton after the splitting. 
What is needed in the NLO-version of CKKW-L in this case is the integral for a 
given factorization scale, which is obtained by simply changing the a s -weight in 
eq. (HU to 

a B (/J.R) fb{x/z,n F ) f a {x,t) 



5 (t) fb{x/z,t) f a (x,H F )' 



(30) 



where z is the energy fraction of the vetoed generated splitting. The derivation in 
eq. (|26[) becomes a bit more cumbersome, but is straight forward. 
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4 Reweighting competing processes 

Often we have many different competing splitting processes. The example we shall 
use here is the process of a quark radiating a gluon {Tg) competing with the process 
of the same quark radiating a photon (-T 7 ). Since generating the latter much less 
likely because of the smallness of qem as compared to a s , the generation may 
become very inefficient if we are interested in observables related to an emitted 
photon. 

In principle we could again consider introducing a boost factor C > 1 and 
replace -T 7 (i) with i%(t) = Cr 7 (t) and do the same with the overestimate t 1 . As 
long as P 1 {t) <C r g {t) we can reweight each event containing n photons with a 
factor 1/C 71 and get approximately the correct results for the observables. How- 
ever this only gives the right emission probability, not the correct no-emission 
probability. 

Instead we adopt different procedure. Every time we generate a photon emis- 
sion (accepted with probability P 7 /_T 7 ), we veto it anyway with a probability 0.5. 
If we veto it, we also reweight the whole event with a factor 2 — 2/C, while if we 
keep it, we reweight the whole event with a factor 2/C. Clearly the emissions will 
still be correctly weighted, 0.5 x 2/C, but now we also get the correct no-emission 
probabilitiefl Loosely speaking we are half the time reweighting the event to com- 
pensate for the boosting of the emission, and half the time compensating for the 
corresponding underestimate of the no-emission probability. 

To see this, we again look at all possible ways of not emitting anything between 
two scales, given by the modified no-emission probability 

A(t , t c ) = A g (t , t c )A 1 (t , t c ), (31) 

where 

~ , - f'° r^(t)dt , . 

Ay(*o,tc) = e J«c ^ (32) 

and the product of weights from all intermediate photon emissions which have 
been vetoed (with probability 1/2, assuming we have already taken care of the 
acceptance factor r-y/T-y): 

°° 9 1 / f to 1 

(w)A(t ,t c ) = £>-~) n i(t ,tc)~ f / dt-r^t) 



71 = 



Au " tc) h(l t ,h 

to 



to 



1 



r 7 (t) 



= A(t ,t c )^ (J °dt [r 7 (*)-A(*)]) 

= A(t ,t c ) (33) 

We note that we could, of course have replaced the probability one half with any 
b, vetoing the emissions with probability b and reweighting with (1 — l/C)/6, while 
reweighting with !/((! — b)C) if not vetoed, and still obtain the correct result. 



^^Note that the whole procedure in principle can be implemented in PYTHIA8 in a non-intrusive 
way, by artificially increasing «em an< i implementing the reweighting and extra rejection in a 
UserHooks class. 
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We also note that while solving the same problem as was addressed in [8] , the 
solution is technically different. In that algorithm, only the standard overestimate 
_T 7 is multiplied by a factor C, while the acceptance of a generated emission at 
scale t is still done with probability _T 7 (t)/_T 7 (t) and a rejected emission instead 
reweights the event by a factor 

A(*)-A(*)/g (34) 
r 7 (t)-r T (t) 

The end result is the same as the method presented here. In fact, if one could 
choos^l a P 1 = 2F 1 and let C — > C/2, the reweighting of the events would be 
exactly the same. 

While we gain in efficiency for the emissions, we will also lose in precision for 
the no-emission probability due to fluctuating weights. It is easy to calculate the 
variance in the weights, but it is maybe more instructive to look at a real example. 

As an illustration we let PYTHIA8 generate standard LEP event, with photon 
emission included in the shower, and we compare the default generation with 
weighting the photon emission cross section with a factor C . We show for different 
C the effect of using the full reweighting procedure {proper weighting), but also 
show for comparison the case of just using event weights with a factor 1/C""' 
{naive weighting). 

In figure [T^ we show the transverse momentum distribution of the most en- 
ergetic photon in an event using C = 1 (i.e. the default), C = 2 for the naive, 
and C = 4 for the proper weighting. The error bands indicates the statistical 
error using 10 8 events, and the results are shown as a ratio to the result from a 
high statistics run (3 ■ 10 9 events) with PYTHIA8. We see that the statistical error 
is somewhat reduced in the reweighted samples, but we also see what seems to 
be a systematic shift in the naive reweighting, due to the mistreatment of the 
no-emission probability. This shift becomes very pronounced if we increase C, as 
seen in figure [T]d, where we use C = 32 for the naive and C = 64 for the proper 
reweighting. Here we see that the statistical errors are very much reduced for 
both reweightings, but the naive procedure is basically useless due to the large 
systematic shift. 

If we require two photons in each event, the gain from the reweighting becomes 
more obvious. In figure [2] we show the distribution in invariant mass of the two 
most energetic photons in an event. Here the gain in statistics is significant also 
for the case of modest boost factors (a), and for the large boost factors in (b) the 
gain in statistics is enormous, while, again, the naive reweighting suffers a large 
systematic shift. 

To isolate the effect on the no-emission probability, figure [3] shows the inclusive 
thrust distribution for the same runs as before. Here we see that especially with 
forceful proper reweighting the statistical error is increased because of the fluctu- 
ating weights, and we see again that the naive reweighting will give a systematic 
shift due to the mistreatment of the no-emission probability. 

2 Note that one need to choose a P which is everywhere some factor higher than _T 7 since 
otherwise the denominator in eq. (1 34 It could tend to zero, giving wildly fluctuating weights. 
3 The proper reweighting has a twice as high boost factor, to get the same weighting of the 
emissions. 
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Fig. 1 The transverse momentum (w.r.t. the thrust axis) of the most energetic photon, given 
as a ratio to the result of a high statistics run with default PYTHIA8 result (see text). 
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Fig. 2 The invariant mass of the two most energetic photons, given as a ratio to the result of 
a high statistics run with default PYTHIA8 result (see text). In (b) the ratio is w.r.t. the result 
for the proper C = 64 reweighting as even with 3 • 10 9 events, the statistical error from the 
default PYTHIA8 run is too large in comparison. 
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Fig. 3 The inclusive thrust distribution, given as a ratio to the result of a high statistics run 
with default PYTHIA8 result (see text). 
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So far we have implicitly assumed that C > 1, since we motivated the whole 
procedure by the desire to increase the number of rare splittings. Note, however, 
that the proof of the procedure does not at all depend on the size of C . In fact it 
can even in some cases be taken negative. 

Consider a case where there are negative contributions to the total splitting 
probability. One of the most simple cases is the emission of a second gluon in a 
final-state dipole shower in e + e~-annihilation into jets. Once a gluon has been 
radiated from the original qq pair, it can be shown that the distribution of a 
second gluon is well described by independent emissions from two dipoles, one 
between the quark and the first gluon and one between the gluon and the anti- 
quark. However, examining the e + e - — > qggq matrix element one finds that there 
is a colour-suppressed negative contribution corresponding to emissions from the 
dipole between the q and q. This contribution is normally ignored completely in 
parton showers, mainly because it is difficult to handle in a probabilistic way in 
the SVA. It may even result in a no-emission probability above unity. 

In the reweighting scheme introduced here we can easily include the negative 
contribution to the splitting functions, and apply a boost factor of C = — 1 for 
the qq-dipole. If a gluon emission is generated from such a dipole, it is then either 
accepted and the event is given a negative weight, or it is rejected, in which case the 
event weight is multiplied by a factor four. We note that in this way it is in principle 
conceivable to implement a parton shower which includes all possible interference 
effects. We will, of course, have even larger issues with statistics, compared to the 
photon emission case above, as we now have potentially large weights that must 
cancel each other, but this procedure could still be an interesting alternative to 
the ones presented in [3] and [J ( an extension of [8] to negative weights) . 

5 Conclusions 

This article does not claim to present innovative new physics results. Rather it 
presents a number of methods collected by the author during a couple of decades 
working with parton showers in general and with the Sudakov veto algorithm in 
particular. They are presented here in the hope that they may come in handy 
for the community now that more and more efforts are put into the merging and 
matching parton showers with matrix element. Especially in the case of matching 
with next-to-leading order matrix elements (and beyond), a thorough understand- 
ing of how parton showers work and knowledge of how to manipulate them is 
necessary, and these kinds of methods may become increasingly important. 
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